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Abstract —A maximum likelihood estimator for fusing the 
measurements in an inertial sensor array is presented. The 
maximum likelihood estimator is concentrated and an iterative 
solution method is presented for the resulting low-dimensional 
optimization problem. The Cramer-Rao hound for the corre¬ 
sponding measurement fusion problem is derived and used to 
assess the performance of the proposed method, as well as to 
analyze how the geometry of the array and sensor errors affect 
the accuracy of the measurement fusion. The angular velocity 
information gained from the accelerometers in the array is shown 
to be proportional to the square of the array dimension and to 
the square of the angular speed. In our simulations the proposed 
fusion method attains the Cramer-Rao bound and outperforms 
the current state-of-the-art method for measurement fusion in 
accelerometer arrays. Further, in contrast to the state-of-the- 
art method that requires a 3D array to work, the proposed 
method also works for 2D arrays. The theoretical findings are 
compared to results from real-world experiments with an in- 
house developed array that consists of 192 sensing elements. 

I. Introduction 

M otion sensing is an essential capability in many sys¬ 
tems to achieve a high level of autonomy. Nowadays, 
inertial motion sensors are ubiquitous in everything from 
industrial manufacturing equipment to consumer electronic 
devices. This widespread usage of inertial sensors has be¬ 
come possible thanks to the last decade’s micro-electrical- 
mechanical-system technology development, which has rev¬ 
olutionized the inertial sensor industry HI. Today, inertial 
sensors can be manufactured at unprecedented volumes and 
at low prices El; six degrees-of-freedom inertial sensor as¬ 
semblies where sold to large-volume customers at less than 
a dollar in 2013 El- Even though development has pushed 
the performance boundaries of the inertial sensor technology, 
contemporary ultralow-cost sensors still cannot fully meet 
the needs of many applications; especially applications where 
the sensor signals have to be integrated over time. These 
applications still suffer from the bias instability, nonlinear¬ 
ities, and thermal instability of the contemporary ultralow- 
cost sensors HI. However, by capitalizing on the decreasing 
price, size, and power-consumption of the sensors, it is now 
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Fig. 1. An in-house designed embedded system with an inertial sensor 
array, which is available under an open-source licence. The array consists 
of 32 MPU9150 Invensense inertial sensor chipsets, each containing an 
accelerometer triad, a gyroscope triad, and a magnetometer triad. 

feasible to construct arrays with hundreds of sensing elements, 
and digitally process these measurements to create virtual 
high-performance sensors. The benefit of using an array of 
sensors is not only an increased measurement accuracy, but 
also an increased operation reliability thanks to the possibility 
of carrying out sensor fault detection and isolation a, Q. 
Further, from the redundant measurements the covariance of 
the measurement errors can be estimated and used to determine 
the reliability of the measurements fbl. Moreover, as will be 
shown, the angular acceleration can be directly estimated and 
for some array configurations the dynamic measurement range 
can be extended beyond that of the individual sensors by utiliz¬ 
ing the spatial separation of the sensors. An example of an in- 
house developed embedded system holding an inertial sensor 
array constructed from contemporary ultralow-cost sensors is 
shown in Fig.[T] The system has 192 inertial sensing elements, 
a micro-controller for information fusion, and a Bluetooth 
interface for communication. Refer to www.openshoe.org and 
E| for details. 

Inertial sensors are primarily used to measure the motion of 
objects. Since a rigid body in a three-dimensional space has 
six degrees-of-freedom, three rotational and three translational, 
the motion of an object is commonly measured using a 
sensor assembly of three gyroscopes and three accelerometers, 
a.k.a. an inertial measurement unit (IMU). However, as a 
rigid object’s motion can also be described in terms of the 

Reproducible research: The layout files and software for 
the inertial sensor array used in the experiments are available 
under an open-source licence at www.openshoe.org. 
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translational motion of three non-collinear points on the object, 
it is also possible to measure the motion using only an assem¬ 
bly of spatially separated accelerometers, a.k.a. a gyroscope- 
free IMU; the object’s rotational motion is then described 
by the relative displacement of the points. Thus, with an 
inertial sensor assembly that consists of both gyroscopes and 
spatially distributed accelerometers, a.k.a. an inertial sensor 
array, rotational information can be obtained both from the 
gyroscopes and the accelerometers. The question is, how 
should the measurements in such sensor array be fused? 

Array signal processing and measurement fusion in the 
context of sensor arrays measuring wave fields is a long- 
studied subject in the signal processing literature la, and 
the within the area developed methods have been applied to 
various types of sensor arrays such as antenna arrays Ha,®, 
magnetic sensor arrays in, acoustic sensor arrays m-m, 
and chemical sensor arrays m- However, since an inertial 
sensor array do not measure a wave field, but something that 
can be referred to as a force field, the signal models and 
methods in the array processing literature are not directly 
applicable to the inertial sensor array measurement fusion 
problem. Nevertheless, based on results from classical me¬ 
chanics about forces in rotating coordinate frames a signal 
model for the inertial sensor array measurements can be 
formulated, and results from the field of estimation theory 
applied to derive an measurement fusion method. Accordingly, 
we will propose a maximum likelihood approach for fusing the 
measurements in an inertial sensor array consisting of multiple 
accelerometer and gyroscope triads. With suitable projection 
vectors, the approach can also be generalized to a system of 
single sensitivity axis sensor units. 

A. State-of-the-art techniques 

Next, a brief summary of existing work on inertial sensors 
arrays is presented. For an in-depth review of the topic the 
reader is referred to the literature survey presented in [Th). 

In the 1960s, it was shown that the specific force and angular 
acceleration of an object can be estimated with an assembly 
of nine accelerometer^ ifTOl : the angular velocity can only 
be estimated up to a sign ambiguity since the centrifugal 
force depends quadratically on the angular speed Il20l - ll22l . 
Though the motion of an object can be tracked using only 
the estimated specific force and angular acceleration, the 
increased sensitivity to measurement errors caused by the 
extra integration needed to calculate the orientation of the 
tracked object renders it as an intractable solution. Therefore, a 
variety of methods to resolve the sign ambiguity of the angular 
velocity estimates by tracking the time development of the 
estimated quantities have been proposed 1201 , l23l - l29l . 

Instead of directly estimating the angular velocity and 
angular acceleration from the accelerometer measurements, it 
is common to first estimate the angular acceleration tensor (or 
variations thereof), and then calculate the angular velocity and 

'with six accelerometers the angular velocity cannot be estimated from 
the measurements at a single time instant, however, a differential equation 
describing the rotation dynamics can be posed by which the angular velocity 
can be tracked over time, see e.g. El, HI 


angular acceleration lag, 1301. To simplify the estimation of 
the angular acceleration tensor, the fact that it only has six- 
degrees of freedom is generally neglected, so that it can be 
estimated via a linear transformation of the accelerometer mea¬ 
surements, see e.g. ifTTl . ial-laa, EI-EI or the Appendix. 
The advantages of this approach are its simplicity and that it 
is straight-forward, using the least squares method (see e.g. 
1281 ). to include the measurements from additional redundant 
accelerometers. The disadvantages are mainly caused by the 
neglected constraints on the tensor, these are as follows: (i) a 
minimum of twelve, instead of nine, accelerometers are needed 
for the estimation problem to be well defined l30l : (ii) the 
estimation accuracy is reduced; and (iii) the locations of the 
sensors must span a 3D space instead of a 2D space l25l . 
The requirement that the locations of the sensors must span a 
3D space significantly increases the size of such a system and 
causes a problem if a sensor array is to be constructed on a 
printed circuit board. 

Since the angular velocity of all points on a rigid object 
are the same, no additional information, except that from 
redundant measurements, is obtained from having several 
spatially distributed gyroscopes in an inertial sensor array. 
Still, a non-negligible reduction in the measurement errors 
can be obtained from the redundant measurements, see e.g. 
the Allan variance analysis in Q. A few methods to fuse the 
measurements from multiple gyroscopes using different filter 
frameworks and signal models are described in E41 - E61 . 

With respect to inertial arrays that consist of multiple IMUs, 
the measurement fusion problem has mainly been studied 
in the framework of global navigation satellite system aided 
inertial navigation systems ©, ED-llQl. In the literature, 
the proposed information fusion approaches can be broadly 
grouped into two categories. In the first category, see e.g. 
0, isi-Ea, the IMU measurements are fused before they 
are used in the inertial navigation systems; commonly a 
weighted average is used and the spatial separation of the 
sensors is neglected. In the second category, see e.g. EOl, 
they are fused after being processed by several parallel inertial 
navigation systems. Refer to 123 and the references therein 
for a discussion on the pros and cons of the two approaches, as 
well as an evaluation of different measurement fusion methods, 
including a few where the spatial separation of the sensors are 
taken into account. 

B. Presented work and findings 

In this paper, the problem of fusing the measurements from 
an array of accelerometer and gyroscope triads is posed as a 
parameter estimation problem for the first time. Results from 
the field of estimation theory are used to derive a maximum 
likelihood based measurement fusion method. The Cramer- 
Rao bound (CRB) for the corresponding estimation problem 
is also derived and used to evaluate the performance of the 
proposed measurement fusion method. Simulations show that 
the proposed method attains the CRB and outperforms the 
current state-of-the-art method for measurement fusion in 
accelerometer arrays. Further, by studying the properties of 
the signal model it is shown that necessary and sufficient 
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conditions for the measurement fusion problem to be well 
conditioned is that the array consists of at least one gyroscope 
triad and three accelerometer triads whose locations are non- 
collinear; current state-of-the-art information fusion methods 
require 3D arrays. Moreover, an analysis of the CRB reveals 
that: (i) the angular velocity information gained from the 
accelerometers is proportional to the square of the array 
dimensions and to the square of the angular speed; (ii) the 
accuracy with which the angular acceleration can be estimated 
decreases with the angular velocity; and (iii) there exists no 
measurement fusion method that can calculate an unbiased, 
finite variance, estimate of the angular velocity for small 
angular velocities from an array of only accelerometer. To 
support the theoretical findings, the accuracy of the proposed 
measurement method is also experimentally evaluated. The 
experimental results show that proposed method works, but 
it also shows that there is a discrepancy between the theo¬ 
retical and experimental results. The discrepancy between the 
theoretical and experimental results are likely primarily due to 
uncertainties in the location of the sensing elements of the real- 
world array. To summarize, the proposed information fusion 
method outperforms the current state-of-the-art method and 
also works for 2D arrays, allowing for smaller and cheaper 
senor arrays to be constructed. However, more research on 
sensor array calibration is needed before the inertial sensor 
arrays can reach their full potential in practice. 

II. Measurement eusion 

In this section a maximum likelihood estimator to fuse 
the measurement from an inertial sensor array consisting of 
Ns accelerometer triads and gyroscope triads will be 
presented. The derivation of the estimator takes its starting 
point in the, from classical mechanics obtained, relationship 
between forces in rotating coordinate frames. 

A. Accelerometer signal model 

Basic kinematics dictates that the specific force in one 
point of a rotating coordinate frame may be decomposed into 
that of another point, a centrifugal term, and an Euler term. 
Specifically, as illustrated in Fig. ID the specific force s^, sensed 
by the uth accelerometer triad located at position in the 
array coordinate frame and which sensitivity axis are aligned 
with the coordinate axis of the array framqfl is given by ll42l 
p. 90] 

Si=s-|-a;x(a;xri)-|-d;xri. (1) 
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Fig. 2. An illustration of the decomposition of the force sensed by an 
accelerometers in an inertial sensor array, overlaid on a picture of the printed 
circuit board of the inertial sensor array shown in Fig. ^ The specific force 
sensed by the fith accelerometer triad is the sum of the specific force at the 
origin of the array coordinate frame, the centrifugal force, and the Euler force. 


Next, let rJa denote the skew-symmetric matrix representa¬ 
tions of the three element vector a for which Jlab = a x b. 
Then ([T]) may be rewritten in terms of matrix multiplications 
as 

= s -I- 

where in the second equality it was used that axb = —bxa. 
The measurement vector y^, consisting of the measurements 
from the accelerometers in the array, can thus be modeled as 

Ys = hs(a;)-f rig, (3) 


where 



1 

e 

I-S 

_1 



-fir/ 

hs(u;) = 


, Hs = [G 0 I 3 ] 

G = 








and cf) = Here and I„ denote a column 

vector of length n with all of the entries equal to one and 
an identity matrix of size n, respectively. Further, rig denotes 
the measurement error of the accelerometers. 


B. Gyroscope signal model 

The angular velocity sensed by the gyroscopes is inde¬ 
pendent of their location in the array. Consequently, the 
measurement vector y^j, consisting of the measurements from 
all of the gyroscopes in the array, can be modeled as 

Yo; = h<^(a;)-f n„, (4) 


Here s denotes the specific force sensed at the origin of 
the array frame. Further, u> and lo denote the array frame’s 
angular velocity and angular acceleration w.r.t. inertial space, 
respectively. Moreover, a x b denotes the cross product 
between vector a and b. This relates the measurements of 
different accelerometers of known locations within the array 
to a common specific force s via the angular velocity a; and 
acceleration uj of the array. 

calibration method to align the coordinate axis of the sensors in an 
inertial sensor anay can he found in ET] 


where h„(a;) = (8 tc’, and denotes the measurement 

error of the gyroscopes. That is, all gyroscopes measures the 
same angular velocity and are disturbed by additive measure¬ 
ment error. 

C. Array signal model 

Concatenating all of the measurements into a single vector, 
the signal model for the full array becomes 

y = h(a;) -f H ^ -f n, 


(5) 
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where 



h(w) = 

hs(tu) 

H = 

H, 

_03Ar^,6_ 

n = 

ns 


Here 0„ „j denotes a zero matrix of size n by m. Thus, the 
signal model (|5]l for the inertial sensor array consists of a non¬ 
linear part h(a;), depending only on the angular velocity us, 
and a linear part H (p, depending only on the specific force s 
and angular acceleration uj. This separation of dependencies 
will turn out to be useful in the estimation of the specific force 
and the angular velocity. 


D. Identifiability conditions 

When = 0, i.e., the array consists of only accelerome¬ 
ters, then h(a;) = hs(a;) = hs(—u;) and the parameters in the 
signal model is not identifiable. Further, H will only have full 
rank if the array has at least three accelerometer triads whose 
locations span a 2D space. To see this, assume without loss 
of generality that the origin of the array coordinate system is 
defined so that ri = O 3 i, then 


03,3 I3 

— rJra I 3 


ftr2 


rank(H) = rank( 


) = 3 -I- rank( 


I 3 

,03AL,3 03 Ar ^_3 


^ri\ 


)• 


Next, note that for ^ 0, rank(r2rj = 2 and the null space 
vector of Cln is r^. Thus, only if there are two accelerometer 
triads whose locations (r^ ^ 0 ) and (r^ 7 ^ 0 ) are such 
that Ti ^ P Tj V/3 £ R, then rank(H) = 6 . That is, H has 
full rank if > 3 and rank([ri ... rTv^]) > 2. Since 
the angular velocity is directly measured by the gyroscopes 
this implies that necessary and sufficient conditions for the 
parameters in the signal model to be identifiable is that: (i) 
the array has at least one gyroscope triad, and (ii) at least 
three accelerometer triads whose locations span a 2D space. 


E. Maximum likelihood measurement fusion 

Assuming the measurement error n to be zero-mean and 
Gaussian distributed with the known covariance matrix Q, the 
log-likelihood function for the signal model Q is given by 

=-i||y - h(a;) - Hc/)||^-i-f c, ( 6 ) 

where c is a constant and ||a||^ = a^Ma. The maximum 
likelihood estimator for {u>, cp} is thus given by 

{tu, cj)} = argmax(L(a;, </>)) 

(jj.cb 

= argmin [\\y - h(a;) - . 

Since the signal model (|5]l is partially linear, we may concen¬ 
trate the log-likelihood function by first fixing the parameters 
ui and maximizing the likelihood function w.r.t. the linear 
parameters cp, and then substitute the result back into the 
likelihood function ll43l . For a fix value u)*, the solution to 


O is given by the weighted least squares solution Il44l . That 
is, 

= argmax(L(u)*, </))) 

<A (8) 

= (HTQ-iH)-iHTQ-i(y - h{ij*)). 

Substituting (l8]l back into (|7]) yields the concentrated maxi¬ 
mum likelihood estimator 

UJ = argmax (Lc(‘^)), (9) 

CJ 

where 

Lc(w) = L{uj, p{uj)) = -^||y - h(a;)||p -f c (10) 

and 

P = Q-i-Q iH(H^Q ^H) 1 . ( 11 ) 

Thus, maximizing the concentrated likelihood problem is the 
equivalent to solving a nonlinear least squares problem. Once 
the angular velocity estimate u: has been found, then the 
maximum likelihood estimate cp can be calculated via ®. 

As the error characteristics of ultralow-cost inertial sensors 
generally vary with motion dynamics and environmental fac¬ 
tors such as temperature P31 , the assumption that the error 
covariance matrix Q is known may not be realistic in all 
situations. The covariance matrix may then be parameterized 
and included in the likelihood function, see e.g., Il46l . A few 
other approaches for estimating the measurement error statis¬ 
tics in inertial sensors arrays are discussed in 0 . However, 
throughout the paper we will assume Q to be known with 
sufficient accuracy. 

Note that the proposed estimator, i.e., (l8]l-([II]l, can also be 
derived without the assumption of the Gaussian distributed 
measurement error using the least squares framework. The 
assumption was introduced to provide stringency with the 
assumptions used to derive the CRB presented in Section HVl 

F. Maximizing the concentrated likelihood function 

The concentrated maximum likelihood estimation problem 
in (|9]l may be solved numerically using the Gauss-Newton 
algorithm Il44l . That is, uj can be iteratively calculated as 

Wfc+l =U>fc + (J^PJ;,)-ljTp(y-h(Wfe)), ( 12 ) 

where k denotes the iteration index and Wq is an initial 
estimate of the angular velocity. A good initial estimate ujq 
to the Gauss-Newton algorithm is given by the weighted 
least squares estimates of the angular velocity calculated from 
the gyroscope readings, i.e., ujq = (( 1 ^^, ® I 3 )Qw^( 1 al ® 
13 ))“^( 1 ^ 0 l 3 )Q~^y(^, where denotes the covariance 
matrix of the gyroscopes measurement error n^. The Jacobian 
3h of h(a;) is given by 

J/t = [A(a;,ri)^ ... A(a;,rwJ^ 1 ^^ 013 ]^, (13) 

where A(u, v) = + J^vxu- 

The performance of the proposed measurement fusion 
method is evaluated in Section |V] where it is compared to 
the CRB derived in Section |IV] Next, we will describe how 
the dynamic range of the angular velocity measurements of 
the array can be extended beyond the dynamic range of the 
gyroscopes used in the array. 
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III. Angular Velocity Measurement Range 
Extension 

The measurement range of contemporary ultralow-cost gy¬ 
roscopes is generally limited to some thousands of degrees per 
second, whereas accelerometers with measurement ranges of 
thousands of meters per second square are widely available. 
Thus, when designing sensor arrays to be used in high dynamic 
applications such as biomedical tracking systems for sport 
activities where the angular velocity may exceed thousands of 
degrees per second, see e.g., ISl, HD, it may be difficult to 
find gyroscopes with a sufficient dynamic range. Consequently, 
the gyroscopes in the array sometimes get saturated. However, 
the proposed measurement fusion method can still be used to 
estimate the angular velocity by simply removing the saturated 
gyroscope readings from the measurement model and adapting 
the dimensions of the measurement error covariance matrix 
accordingljH 

When all of the gyroscopes are saturated, then the mea¬ 
surement model Q becomes the equivalent to the model in 
OJ, for which the angular velocity is not identifiable since 
hs(a;) = hs(—tu). However, since the sign of the elements 
in the angular velocity vector can be deduced from the 
saturated gyroscopes, the angular velocity is still identifiable. 
In Fig. [3] an example of the angular velocity measurement 
range extension is shown, when applied to data collected by 
the array in Fig. [T] 

When the gyroscopes are saturated, finding a good initial¬ 
ization of the Gauss-Newton algorithm becomes trickier and 
the convergence of the algorithm can be slow. In applications 
where the inertial sensor array is used on a regular basis 
to measure the specific force and angular velocity, an initial 
estimate may be found by predicting the angular velocity from 
the at pervious time instant estimated angular velocity and 
angular acceleration ED. Another method to find an initial 
estimate when the gyroscopes are saturated and a 3D array is 
used is through the angular acceleration tensor method that is 
summarized in the Appendix. 

The performance of the proposed measurement fusion 
method, for the case when the gyroscopes are saturated, is 
evaluated in Section |V] The performance is compared to the 
CRB derived in the next section. 

IV. Cramer-Rao Bound 

Since the measurement error in the signal model is as¬ 
sumed to be Gaussian distributed, a general expression for 
the CRB for the measurement fusion problem at hand is 
straight forward to derive starting from the results in ll44l p.49]. 
Assuming the following parameter order 6 = [o;^ 0^] = 

[o;^ s^] , the CRB becomes 

Cov(6>) h 2:(0)-\ (14) 

and the Fisher information matrix I{6) is given by 

x{e) = (15) 

^Theoretically, simply removing the measurements that ai‘e believed to be 
saturated will create a bias in the estimated quantities, however our practical 
experience shows that this effect is negligible. 


Accelerometer measurements versus time. 



(a) Overlaid measurements form all of the accelerometers in the array. The 
differences between the measurements are caused by the spatial separation 
of the sensors and can be used to estimate the angular velocity and angular 
acceleration of the aiTay. 


Gyroscope measurements versus time. 



(b) Overlaid measurements form all of the gyroscopes in the array. The 
gyroscopes saturates at 2000 °/s. 


Estimated angulai' velocity versus time. 



(c) Angular velocity estimated using the proposed method. The angular 
velocity can be estimated even though the gyroscopes are saturated. 


Fig. 3. An illustration of the angular velocity measurement range extension 
using the proposed method. The measurements were recorded while twisting 
the array in Fig. ^ by hand. The differences between the accelerometer 
measurements, see Fig. |3(a)| are used to estimate the angular velocity, see 
Fig. |3(c)| even when the gyroscopes are saturated, see Fig. |3(b)| 

where €> = [J/^ H]. Here the notation B V A implies 

that (B — A) is a positive semi-definite matrix. Since the 
measurement model is linear in 0, the Fisher information 
matrix will not depend on the specific force and the angular 
acceleration. To gain a deeper understanding of the estimation 
problem at hand, we will study the CRB for a special set of 
array configurations next. 

A. Square arrays with uncorrelated measurement errors 

Consider the case where the following three conditions 
apply; (i) the accelerometers are mounted in a planar square 
grid with a spacing a in between each sensor and the origin 
of the array coordinate system is defined at the center of 
the grid; (ii) the measurement errors are uncorrelated; and 
(iii) all of the accelerometers and gyroscopes have the same 
error variance. Then Q = (fT^IsArJ © (o’JIsal), where © 
denotes the matrix direct sum. Further, and crj denote 
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the measurement error variance of the accelerometers and 
gyroscopes, res^ctively. Moreover, thanks to the symmetry 
of the array The Fisher information matrix 

then takes the following form 


where 


i{e) 


■Xn 

X12 

03,3 

X12 

X22 

03,3 

03,3 

03,3 

Ns-f 


n- _^T J 1 

Xn - ^l3 + 

UJ ^ s 

Xi 2 = 

X22 = -^G^G. 


(16) 


(17) 

(18) 
(19) 


Here denotes the Jacobian of hs(a;). From the Fisher 
information matrix in ( fTSl l we can see that for the considered 
type of array, the CRB for the specific force is given by 
Cov(s) ^ Xg , where X^ = I3. That is, the accuracy 

with which the specific force can be estimated is inversely 
proportional to the number of accelerometer triads in the array, 
and independent of the angular velocity and geometry of the 
array. This is quite intuitively since for a symmetric array with 
a set of accelerometer triads having the same measurement 
error characteristics (measurement error covariance), the max¬ 
imum likelihood estimator for the specific force is given by 
the arithmetic mean of the accelerometer measurements. 

Using the Schur complement, the CRB for the angular 
velocity can be found to be 


Coy{<:>) h X-i (20) 


where 


and 


X^ —1^11 -^12-^22 -^12 


Ns 

rii = ^ A(a;, r,)TA(w, r,) 

Ns 

ri2=J^^G = -^A(u;,r,)Tnr, 

i^l 

Ns 

r22=G^G = Y,^JNri 

2=1 


( 21 ) 


( 22 ) 

(23) 

(24) 


Next, using the symmetric properties of the considered array 
geometry and performing, tedious but straight forward, calcu¬ 
lations give that 


X — 

-L-uj — 


-f 


6 cr? 


+ Wy 


2ujy + 
2WyUlz 


2uixUJ^ 
2UJyUJ z 


( 25 ) 


Note that the assumption about the accelerometers be¬ 
ing mounted in a planar square grid implies that Ns S 
{4, 9,16, 25,...}. From the expression for the Fisher informa¬ 
tion matrix in (l25T l. we can see that when the array is stationary 
no rotational information is gained from the accelerometers, 
and the CRB becomes equivalent to the covariance of the 
arithmetic mean of the gyroscope measurement errors. When 
the array starts rotating, the accelerometers will start to pro¬ 
vide rotational information and the accuracy with which the 
angular velocity can be estimated increases proportionally to 
the squared angular rates. Further, the rotational information 
gained from the accelerometers scale quadratically with the 
distance a between the sensors. 

When all gyroscopes are saturated, i.e., \uji\ > 7 , V* S 
{x,y^z}, where (— 7 , 7 ) is the dynamic range of the gyro¬ 
scopes, no information other than the sign of the angular 
velocity is provided by the gyroscopes and the angular ve¬ 
locity must be estimated sole from the accelerometers. As the 
sign information can be seen as a special case of inequality 
constraints, and inequality constraints generally does not have 
any effect on the CRB B9l . the Fisher information for this case 
can be found by letting crj — 00 . The inverse of the Fisher 
information matrix (l25l l then takes the following simple form 


1 (|.-._|>7) 


a^{N^-Ns) 


0 






(26) 


where the non-zero of the diagonal elements have been left 
out for brevity. From this, we can see that in an array where 
the gyroscopes have an infinitely small dynamic range, i.e., the 
gyroscopes only provides information regarding the sign of the 
angular velocity, then the covariance of the angular velocity 
estimates tends toward infinity as the angular velocity goes 
toward zero. Hence, for so called gyroscope-free IMUs, i.e., 
an array of only accelerometers, there is no estimator that can 
provide an unbiased, finite variance, estimate of small angular 
velocities. Therefore, the practical use of gyroscope-free IMUs 
in inertial navigation systems for low-dynamical applications 
is questionable. 

The CRB for the angular acceleration 6j is given by 


Cov(cI>) hXvi, (27) 
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From (|28 T i we can note that the accuracy with which the 
angular acceleration uj can be estimated, somewhat surpris¬ 
ingly, decreases with an increasing angular velocity. The best 
accuracy is achieved when the array is in linear motion 
(uj = 0 ) and is given by 
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Fig. 4. Geometries of the accelerometer triads used in the Monte Carlo 
simulations. 

V. Simulation and Experiments 

In this section the performance of the proposed maxi¬ 
mum likelihood measurement fusion method will be evaluated 
through simulations and real-world experiments. The results 
will be compared to the derived CRB. 

A. Simulations 

To evaluate the accuracy of the estimator we have con¬ 
ducted two Monte-Carlo simulations, where each Monte-Carlo 
simulation runs the proposed measurement fusion method on 
10^ data realizations. The array considered in the simulations 
consists of four accelerometer triads and four gyroscope triads, 
where the accelerometer triads are mounted as illustrated in 
Fig. |4(a)| Note that all of the arrays that fulhll the assumptions 
made in Section IIV-AI can be transformed into an equivalent 
array with a geometry as in Fig. |4(a)| plus possible additional 
sensors located at the origin. The measurement errors of the 
accelerometers and the gyroscopes were assumed to be uncor¬ 
related and to have the standard deviation a1 = 0.01 [m/s^] 
and CTtj = 1 respectively. Further, the gyroscopes were 

assumed to saturate at 2000 [°/s^]. These parameter values 
were selected to reflect the performance that can be expected 
from ultralow-cost inertial sensors during high dynamic opera¬ 
tions when scale-factor, g-sensitivity, and cross-axis sensitivity 
errors become the main error sources. The small value selected 
for the accelerometer measurement error variance is motivated 
by the fact that the accelerometers in the array can easily 
be calibrated using e.g., the method in HTI . A calibration of 
the gyroscopes is complicated and requires a rotational rig. 
If the accelerometers where uncalibrated, then the standard 
deviation of the accelerometer measurement errors would be 
a magnitude higher approximately; refer to ED for typical 
error figures before and after sensor calibrations. 

The root mean square errors (RMSE) calculated from the 
Monte-Carlo simulations, together with the corresponding 
square root of the CRBs, are shown in Fig. m In Fig. 1^ 
the result when the array is rotated around the x-axis, i.e., 
an in-plan rotation, is shown. Clearly, no information about 
the angular rate around the z-axis (pointing out of the plan) 
is gained from the accelerometers, whereas the accelerometer 
measurements contribute to the estimation of the angular rates 
around the x-axis and y-axis. In Fig. |5(b)| the result when 
the array is rotated around the z-axis, i.e., an out-of-plan 
rotation, is shown. In this case no information about the 




Eig. 5. The angular velocity estimation RMSE of the proposed measurement 
fusion method for a planar array with the geometry illustrated in Fig. |4(a)| 
The estimation accuracy obtained by simply averaging the gyroscope mea¬ 
surements is given by the gray horizontal line. 


rotations around the x-axis and y-axis are gained from the 
accelerometers; only information about the rotation around the 
z-axis is gained. The estimation accuracy achieved by simply 
averaging the gyroscope measurements is illustrated by the 
grey horizontal line in the hgures. 

To compare the proposed measurement fusion method with 
the in gyroscope-free IMUs commonly used angular accelera¬ 
tion tensor method, a Monte Carlo simulation using the array 
geometry illustrated in Fig. |4(b)| was also conducted; the angu¬ 
lar acceleration tensor method is summarized in the Appendix. 
The non-planar array geometry is needed for the tensor method 
to work. The result of the simulation is shown in Fig. |6] from 
which it can be seen that the tensor method does not achieve 
the CRB and is outperformed by the proposed method. Note 

















































Maximum likelihood versus tensor method. 
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Fig. 6. The angular velocity estimation RMSE of the proposed maximum 
likelihood measurement fusion method and the angular accelerati on ten sor 
method. The geometry of the considered array is illustrated in Fig. |4(h)| and 
the array is rotating around all of the axes, i.e., lj = ||aj||/\/3 [ 111 ]^. 
The estimation accuracy obtained hy simply averaging the gyroscope mea¬ 
surements is given by the gray horizontal line. 



Fig. 7. RMSE of the from real-world data calculated angular speed. Also 
shown, the RMSE of the estimated angular speed when calculated from 
simulated data that includes random sensor location errors. 


revision of the array when used in a foot-mounted inertial 
navigation system were evaluated. 


that the tensor method only uses the measurements from the 
accelerometers and the sign information from the gyroscopes, 
and therefore the behavior of the two methods should only 
be compared for the part of the simulation where all of the 
gyroscopes are saturated, i.e., above 2000 o /s. 

B. Experiments 

To evaluate the proposed method with real-world data the 
in-house designed inertial sensor array, shown in Fig. [T] 
was mounted in a mechanical rotation rig and data was 
recorded at different angular speeds. At each angular speed, 
data corresponding to 100 time instants were recorded. The 
data was processed using the proposed maximum likelihood 
method with the same noise variance settings as used in the 
simulations and the RMSE of the estimated angular speed was 
calculated; the angular speed and not the angular velocity as 
an evaluation criterion was chosen because of the practical 
problem of accurately aligning the coordinate axes of the array 
with the rotation axis of the rotation rig. The result is displayed 
in Fig. |7l and it shows that the proposed method works but 
it does not achieve the accuracy predicted by the theoretical 
model. Several factors such as sensor scale factor errors, 
misalignment between the sample instances of the sensors, 
and uncertainties in the location of the sensing elements, 
have been identihed as likely sources for the discrepancy. For 
illustrational purposes the result of a Monte-Carlo simulation 
resembling the experimental setup, and where random errors 
with a standard deviation of 0.1 mm were added to the sensor 
locations, is also shown in Fig. |70 More experimental results 
can be found in ll50ll . where the performance of an earlier 

^Placement errors in the order of 0.1 mm are likely to arise in the PCB 
manufacturing and population process. 


VI. Summary, Conclusions, and Future Research 

Approximately 5 billion smartphones were sold in the 
hve year span 2011-2015, of which most where equipped 
with some type of inertial motion sensors. Thus, the smart¬ 
phone industry has become a driving force in the devel¬ 
opment of ultralow-cost inertial sensors, and six-degree-of- 
freedom inertial sensor assemblies are sold to large-volume 
costumers for less than a dollar. Unfortunately, these ultralow- 
cost sensors do not yet meet the needs of more demanding 
applications like inertial navigation and biomedical motion 
tracking systems. However, by adapting a wisdom of the 
crowd’s thinking and design arrays consisting of hundreds of 
sensing elements, one can capitalize on the decreasing cost, 
size, and power-consumption of the sensors to construct virtual 
high-performance low-cost inertial sensors. Accordingly, we 
have proposed a maximum-likelihood method to fuse the 
measurements in arrays consisting of multiple accelerometer 
and gyroscope triads. The proposed method has been eval¬ 
uated through simulations as well as real-world tests with 
an in-house designed array with 192 sensing elements. The 
simulation results show that the proposed method attains the 
CRB and it outperforms the current state-of-the-art method. 
Further, compared to the state-of-the-art method the proposed 
method also works with 2D arrays, which is fundamental for 
the production of arrays on a single printed circuit board. 
Moreover, by utilizing the spatial separation between the 
accelerometers, the angular velocity measurement range can be 
extended beyond that of the gyroscopes. The experimental re¬ 
sults show that the proposed methods work, but do not achieve 
the accuracy predicted by the theoretical model. Uncertainties 
in the position of the sensing elements have been identihed as a 
likely source for the discrepancy. Further research is needed to 
verify the source of the problem, and new inertial sensor array 
calibration methods must be developed. Moreover, information 
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fusion methods that also considers the time development of 
the sensor signals and allows for simultaneous calibration 
and motion tracking should be investigated. Another open 
research question is to design a measurement error covariance 
estimation method that enables automatic weighting of the 
sensor measurements. 

Appendix 

The angular acceleration tensor method used in the simu¬ 
lations are here summarized. For a more detailed description 
see ll28ll . Introducing the angular acceleration tensor W = 
and using the properties of the cross product, the 
model for the specihc force (|2]i can be rewritten as 




— Hr 


Wr,. 


(30) 


The output of the Ng accelerometer triads can be described 
by the following linear matrix equation 


where 


Y = [ y., . 
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■ 1 ... 
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Y = XR + N (31) 

(32) 

( 33 ) 

Here and ris^ denote the measurement and measurement 
error of the i:th accelerometer triad. Neglecting the fact that 
the tensor W only has six degrees of freedom, the least square 
estimate of the matrix X is given by 


X = YRT(RR 


T\-l 


(34) 


Noteworthy is that since X is a 3-by-4 matrix, the mea¬ 
surements from at least four accelerometer triads are needed. 
Further, for R to have a full row rank and the estimation 
problem to be well defined, the locations of the accelerometer 
triads must span a three dimensional space. From the estimated 
tensor W, the angular acceleration can be calculated as 




^ = [ t«3,2 - W2, 


3 - ^3,1 W2,l - tUl,2 


(35) 


where Wi^j = The angular velocity can, up to a sign 

ambiguity, be calculated from the left hand side of the equality 


= i(W- 


W 


)-itr(W- 


'W^)l3 


(36) 


Here tr( ) denotes the trace operator. 
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